Su un'estensione della teoria di Lagrange 
per i moti secolari. 

Antonio Giorgilli*, Ugo Locateli^, Marco Sansottera* 



Sommario 

La teoria di Lagrange per i moti secolari delle eccentricità ed in- 
clinazioni delle orbite planetarie si fondava su un'approssimazione, 
dettata in larga misura dalla complessità dei calcoli necessari, che 
consisteva nel considerare solo equazioni lineari. In questa memoria 
riprendiamo in considerazione i metodi di Lagrange alla luce della teo- 
ria della stabilità esponenziale di Nekhoroshev. Grazie agli algoritmi 
sviluppati negli ultimi anni e alle tecniche di manipolazione algebrica 
possiamo tener conto anche dei contributi non lineari alle equazioni. 
Come applicazione cerchiamo di determinare i tempi di stabilità per 
il problema dei tre corpi nel caso del Sole e dei due pianeti maggiori, 
Giove e Saturno, mostrando che si possono ottenere risultati realistici, 
ancorché non ottimali. 

Lagrange's theory for the secular motion of perihelia and nodes of 
the planetary orbits was based on consideration of a linear approssi- 
mation of the dynamical equations, compatible with the complexity of 
the calculations. We extend Lagrange's investigations in the light of 
Nekhoroshev's theory of exponential stability. Using effective algori- 
thms recently developed and computer algebra we investigate the non 
linear problem. We apply our methods to the problem of three bodies 
in the Sun-Jupiter-Saturn case, thus showing that realistic results, 
although not optimal, can be obtained. 
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1 Introduzione 



Nel 1782, in due corpose memorie presentate all'Accademia di Berlino, Giu- 
seppe Luigi Lagrange pubblica la versione estesa della sua teoria sui mo- 
ti secolari dei nodi e dei perieli dei pianeti, e conclude il suo studio con 
l'affermazione della stabilità del Sistema Solare. Si tratta, per quel tem- 
po, di un risultato di notevole rilevanza, soprattutto se si tiene conto che 
viene pubblicato in un periodo in cui sono molto vive le discussioni sulle 
"ineguaglianze secolari", in primis quella di Giove e Saturno, e sull'effettiva 
possibilità di spiegare tutte le apparenti irregolarità dei moti planetari sulla 
base della gravitazione newtoniana. 

Alla luce delle nostre conoscenze attuali, e con la sensibilità matemati- 
ca del nostro tempo, quel risultato non appare completamente rigoroso, in 
quanto frutto di approssimazioni la cui completa validità non è assicurata. 
In questa memoria vogliamo illustrare come ed in che senso si possa esten- 
dere la teoria di Lagrange tenendo conto degli sviluppi recenti delle nostre 
conoscenze. 

Il nostro metodo si ricollega in modo diretto a quello di Lagrange in quan- 
to prendiamo come riferimento le orbite circolari, studiando poi l'evoluzione 
delle eccentricità e delle inclinazioni come oscillazioni intorno ad un equili- 
brio. L'estensione rispetto al metodo di Lagrange prende avvio dalla teoria 
della stabilità esponenziale, sviluppata in forma teorica da Moser [TJ], Lit- 
tlewood [Ti] [12] e Nekhoroshev [16] [T7]. Al fine di applicare tale teoria, 
dobbiamo procedere allo sviluppo in serie della perturbazione mutua di Gio- 
ve e Saturno, tenendo conto anche dei termini di secondo ordine nelle masse 
e dei termini non lineari nelle eccentricità. Ciò risulta fattibile grazie alla 
disponibilità dei metodi di manipolazione algebrica al calcolatore, che ren- 
dono calcolabili sviluppi in serie praticamente impossibili in passato. Questo 
aspetto viene discusso brevemente nel paragrafo [2, in cui richiamiamo anche 
i punti essenziali del metodo di Lagrange. 

In un secondo tempo facciamo ricorso al metodo della forma normale di 
Birkhoff rienunciato mediante algoritmi di calcolo esplicitamente applicabili 
grazie alla manipolazione algebrica. Questi metodi vengono richiamati nel 
paragrafo [3j Infine mostriamo come si possano ottenere stime di stabilità 
su tempi lunghi sfruttando la forma normale ad ordine finito; questo viene 
discusso nel paragrafo HI 

Il paragrafo E] riporta i risultati dell'applicazione del nostro metodo al 
sistema Sole-Giove-Saturno. Le conclusioni, che discutiamo brevemente nel 
paragrafo [6j non saranno completamente soddisfacenti: si arriva a garantire 
la stabilità solo sull'arco di IO 7 anni, nettamente inferiore rispetto alla durata 
del Sistema Solare ed anche rispetto ai tempi calcolati mediante integrazione 
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diretta delle equazioni. Si tratta comunque di uno tra i migliori risultati 
oggi disponibili sulla base dell'applicazione dei metodi perturbativi ed in 
particolare mostra come il meccanismo della stabilità su tempi esponenziali 
possa ben essere significativo almeno per i pianeti maggiori del nostro sistema. 



2 II problema secolare 

Lo schema di sviluppo perturbativo che utilizziamo può considerarsi come 
una riformulazione in ambito hamiltoniano dello schema di Lagrange. Que- 
sta parte del calcolo risulta alquanto laboriosa, ma omettiamo i dettagli in 
quanto si tratta di un argomento classico che si può trovare ad esempio nelle 
Legons de Mécanique Celeste di Poincaré. Partendo dall' Hamiltoniana del 
problema dei tre corpi si procede in modo classico introducendo le coordi- 
nate eliocentriche, che consentono di eliminare il moto del baricentro, ed 
effettuando la riduzione del momento angolare. In tal modo si eliminano 5 
gradi di libertà, sicché ci si riduce a considerare 4 sole coppie di coordinate 
canoniche. 

Si considerano poi gli elementi orbitali, ossia il semiasse maggiore a, 
l'eccentricità e, l'inclinazione i, l'anomalia media £, l'argomento del perie- 
lio u e la longitudine del nodo Q, ed in termini di questi si introducono le 
variabili di Poincaré definite come 



Aj = fj.jy/Gimo + m^aj , & = y/2A~^l- ^fl^ 



e"* cos Uj , 



Xj = ij + Uj , rjj — — y/2Aji 1 — A /l — e? sin . 



dove m , m b m 2 sono le masse dei tre corpi, \ij = — ° +r ^ con j = 1,2 le 
masse ridotte dei due pianeti e G è la costante di gravitazione. Qui abbiamo 
omesso le variabili che descrivono le inclinazioni ed i nodi perché vengono 
eliminate dalla riduzione del momento angolare. Ricordiamo che gli angoli 
A e le azioni ad essi coniugate A vengono detti variabili veloci, in quanto le 
frequenze corrispondenti sono quelle del moto kepleriano, mentre le £j, rjj, 
che si riferiscono alle eccentricità, sono dette variabili lente o secolari, in 
quanto nell'approssimazione kepleriana esse restano costanti mentre nelle 
approssimazioni successive sono soggette a variazioni molto lente visibili solo 
sull'arco di secoli. 

L' Hamiltoniana in variabili di Poincaré comprende due contributi 

H = H (A) + H 1 (A,\,Z, V ) , (1) 
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che consistono in una parte imperturbata H (A) che descrive il moto ke- 
pleriano (ellittico o circolare) ed una perturbazione i?i(A, A, £, 77) dovuta 
all'interazione mutua tra i pianeti. L'Hamiltoniana così ottenuta può es- 
sere sviluppata in serie di potenze intorno alle orbite kepleriane circolari. A 
tal fine si sceglie il valore A* risolvendo l'equazione 

= n*, j = l,2, 

A=A* J 
5=77=0 

dove il simbolo (•)> indica la media rispetto agli angoli veloci e n* è la frequen- 
za fondamentale del moto medio relativa all'angolo Xj, per una esposizione 
più dettagliata si veda [T3j . Si procede poi ad uno sviluppo in serie di po- 
tenze nelle variabili £, 77 nell'intorno dell'origine, osservando che £ = rj = 
corrisponde ad eccentricità nulla. Definendo A = A* + A' sviluppiamo in serie 
di potenze di A' l'Hamiltoniana (DQ) ed eliminando gli apici per semplificare 
le notazioni la riscriviamo nella forma 

H (A) = J2^ + O(A 2 ) , (2) 

i 

che descrive il moto circolare con frequenze kepleriane u, e 

H 1= J2 9, fc (A,A)£V, (3) 

che è lo sviluppo della perturbazione in serie di potenze di £, 77 con coefficienti 
c^A, A) periodici in A. 

Lo schema generale esposto fin qui non differisce di molto da quello di 
Lagrange, ma nel nostro calcolo introduciamo due differenze significative. La 
prima è che, grazie all'uso della manipolazione algebrica al calcolatore, siamo 
in grado di calcolare esplicitamente anche molti termini di secondo ordine 
nelle masse e termini non lineari nelle variabili lente. La seconda differenza 
è che possiamo garantire un'approssimazione migliore per la dinamica dei 
semiassi maggiori, proprio tenendo conto della perturbazione fino all'ordine 
2 nelle masse. A tal fine riordiniamo lo sviluppo della perturbazione nella 
formula (jHJ) come 

F X (A, A, £, 77) = / (A, £, v) + A(A, A, £, V ) + 0(A 2 ) , (4) 

dove si intende che /o(A, £,77) e /i(A, A, £, 77) sono rispettivamente di grado 
zero e uno nelle azioni veloci A, e con una coppia di trasformazioni canoni- 
che cerchiamo di rimuovere /o ed fi in modo da lasciare tra i contributi che 
dipendono effettivamente da A solo quelli che sono almeno di ordine 2 nelle 
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masse. Questo procedimento si ispira alla costruzione della forma normale di 
Kolmogorov, ed è descritto in dettaglio in [H]. Ne risulta un'Hamiltoniana 
che ha ancora la forma ([3]), ma la perturbazione Hi non contiene più al- 
cun termine effettivamente dipendente da A che sia anche indipendente da o 
lineare in A e di ordine inferiore a 2 nelle masse. 

Infine, riprendendo lo schema di Lagrange, introduciamo il modello seco- 
lare. A tal fine, con un'operazione di media, eliminiamo dall' Hamiltoniana la 
dipendenza dagli angoli veloci A e fissiamo A. Concretamente ciò si ottiene 
eliminando dagli sviluppi in serie di Fourier tutti i termini che contengono le 
variabili A e ponendo A = 0. Ciò corrisponde a fissare la dinamica dei semias- 
si maggiori in modo che sia una piccola variazione quasi periodica dell'orbita 
circolare. In tal modo l'Hamiltoniana risultante dipende solo dalle variabili 
lente £, rj, e si riduce ad un sistema a due gradi di libertà. Il fatto rilevante è 
che lo sviluppo in serie di potenze dell 'Hamiltoniana non contiene contributi 
di grado dispari, ed in particolare neppure termini lineari, sicché si deve stu- 
diare la dinamica di un sistema conservativo in prossimità di un equilibrio. 
Precisamente si ottiene un'Hamiltoniana della forma 

H(£,n) = H (Z,ri) + H 3 (£ ì ri) + ... , (5) 

dove Ho, i?2, • • • sono polinomi omogenei di grado rispettivamente 2, 4, . . .. 

La proprietà che abbiamo appena enunciato è a prima vista sorprenden- 
te, ed in effetti a Lagrange spetta il merito di averla messa in evidenza per 
primo e di aver fondato su di essa le sue ricerche sui moti secolari. Il proce- 
dimento da lui seguito può riformularsi dicendo che si considera la sola parte 
quadratica dell 'Hamiltoniana, sicché si deve studiare un sistema di equazioni 
lineari a coefficienti costanti. Il metodo per risolvere tali equazioni era ben 
noto a Lagrange, dato che egli stesso lo aveva sviluppato in una memoria del 
1763 dandogli sostanzialmente la forma che ancora troviamo nei trattati di 
Analisi Matematica. 

Il suo risultato di stabilità consiste poi nel mostrare che le soluzioni 
scritte come composizione di moti periodici con le frequenze e le ampiez- 
za calcolate per i pianeti restano sempre limitate in un intorno abbastanza 
piccolo dell'origine. La conclusione si fonda sull'assunzione, non dimostra- 
ta ma accettata come perfettamente plausibile, che la dinamica nell'intorno 
dell'equilibrio non venga influenzata in modo rilevante dai contributi non li- 
neari che compaiono nelle equazioni. Proprio questo invece è il punto che dob- 
biamo rimettere in discussione alla luce degli sviluppi delle nostre conoscenze 
dopo la scoperta dei moti caotici da parte di Poincaré. 

Torniamo dunque a considerare l'Hamiltoniana (J5J). Grazie all'analiticità 
di tutte le funzioni e le trasformazione coinvolte, fin qui lo sviluppo risulta es- 
sere convergente in un intorno dell'equilibrio. Inoltre, avendo determinato le 
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frequenze ooj dei moti secolari col metodo di Lagrange, possiamo anche intro- 
durre una trasformazione lineare di coordinate che pone la parte quadratica 
dell' Hamiltoniana nella forma diagonale 

^) = £f(^), (6) 

3 

3 La stabilità esponenziale 

A questo punto ha inizio la parte più rilevante della nostra estensione del 
lavoro di Lagrange, in quanto teniamo conto proprio dei contributi non lineari 
all' Hamiltoniana ([5$ che abbiamo potuto calcolare grazie alla manipolazione 
algebrica con uso del calcolatore. 

Si fa ricorso alla forma normale di Birkhoff nell'intorno dell'equilibrio. Si 
considera un' Hamiltoniana della forma 

H{x, y) = H (x, y) + H x {x, y) + H 2 (x, y) + . . . , (7) 

dove 

j 

e Hi(x, y), H 2 (x, y) . . . sono polinomi omogenei di grado 3, 4, ... , sicché si ha 
una serie di potenze convergente in un intorno dell'origine. L' Hamiltoniana §5$) 
ha questa forma, salvo la particolarità di non contenere termini dispari nello 
sviluppo. L'obiettivo è costruire una trasformazione canonica di coordinate 
prossima all'identità che ponga l'Hamiltoniana nella forma 

y) = H Q {x, y) + Z 1 (x,y) + ... + Z r (x, y) + TZ^ r+1 \x, y) , (8) 

dove Zi(x, y), . . . , Z r (x, y) dipendono solo dalle azioni Ij = (x^j +y|)/2, e 
è un resto non normalizzato di grado almeno r + 3 nelle variabili 
x, y. A tal fine utilizziamo una successione di trasformazioni canoniche ge- 
nerate mediante l'algoritmo della serie di LieQ Precisamente, supponendo di 
aver costruito la forma normale Z^ r ~ 1 ^ fino ad un ordine r — 1, si determina 
una funzione generatrice Xr{%,y) come polinomio omogeneo di grado r + 2 
risolvendo l'equazione 

LhqXt — Z r = Q r , (9) 

dove Q r è la parte omogenea di grado r del resto non ancora normalizzato, e 
Lf — {/, •} è la parentesi di Poisson con la funzione /, ovvero la derivata di 

^er un'esposizione del metodo delle serie di Lie in ambito hamiltoniano si veda ad 
esempio [S]. 
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Lie lungo il campo hamiltoniano generato da /. Si determina poi la nuova 
Hamiltoniana calcolando 

Z^=ex V (L Xr )Z^ , 

dove 

exp(L Xr )=l + L Xr + ±Ll + ±Ll + ... , 

è l'operatore esponenziale della serie di Lie. Lo schema di calcolo è facilmente 
programmabile mediante un manipolatore algebrico, in quanto richiede solo 
il calcolo di somme, prodotti e derivate di polinomi omogenei. 

Se, ignorando per un momento il problema della convergenza del procedi- 
mento di costruzione della forma normale, si immagina di applicare infinite 
volte lo schema appena descritto si ottiene una forma normale 

Z^\x,y) = H (I) + Z X {I) + Z 2 (I) + . . . , 

funzione solo delle azioni Ij = [x 2 - + y|)/2 che risultano essere costanti del 
moto per Z^°°\ Scrivendo le equazioni di Hamilton si ottiene così 

Xj = £lj(I)yj , i/j = -tljWxj , 

dove 

«,(/> -, + fffl + §fm + - • 

sono le frequenze che sono costanti del moto, essendo funzioni solo delle /, 
e in conseguenza della non linearità delle equazioni dipendono dal valore 
iniziale delle azioni /. Le equazioni sono integrabili in modo elementare, in 
quanto le soluzioni sono oscillazioni con frequenze O(J) dipendenti dal dato 
iniziale. Se così fosse potremmo affermare di aver esteso la teoria di Lagrange 
nel senso che abbiamo calcolato dei valori migliori per le frequenze secolari, 
mantenendo poi tutte le sue conclusioni per quanto riguarda la stabilità. 
Inoltre si giustificherebbe la validità della teoria di Lagrange in quanto per 
piccole ampiezze le correzioni non lineari alle frequenze sono piccole. 

Vi sono però due difficoltà. La prima, nota come problema dei picco- 
li divisori, e che la soluzione dell'equazione <Q e possibile solo assumendo 
delle condizioni di non risonanza sulle frequenze u del moto imperturbato. 
Precisamente, dal punto di vista teorico si chiede che la quantità . kjUij 
con coefficienti kj interi si annulli solo se kj = per tutti i j. La seconda 
difficoltà si cela nella falsità dell'ipotesi che la forma normale risulti essere 
convergente. Non è difficile verificare che ciascun passo del procedimento è 
ben definito, nel senso che la funzione Z^ r \ per ogni r finito, risulta essere 
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convergente in un intorno dell'origine, ad esempio una sfera di raggio p r . Le 
stime analitiche però consentono solo di dimostrare che p r è limitato infe- 
riormente da una successione che tende a zero almeno come 1/r. Questo 
non esclude che si possa avere convergenza in casi specifici, ed in effetti si 
possono costruire esempi di Hamiltoniane che ammettono una forma normale 
di Birkhoff convergente. Tuttavia nel lavoro di Siegel [TB] si mostra che la 
divergenza è il caso tipicoH 

Ciò che rende utile lo sviluppo perturbativo nonostante la divergenza è 
il carattere asintotico delle serie. Consideriamo un intorno A p dell'origine 
imponendo la condizione \lj(x, y)\ < p 2 /2 per j = l,...,n. Assumiamo 
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poi che le frequenze u> soddisfino la condizione diofantea 

con 7 > e r > n - 1, essendo n il numero di gradi di libertà, e \k\ = 
Y2j \kj\. Allora con considerazioni teoriche si può verificare che nell'intorno 
considerato si ha 

sup \ll ir+1) \ < B r (r\) n+1 p r+1 , 

dove B è una costante positiva e n è il numero di gradi di libertà del sistema 
(si veda ad esempio j3] o jl]). Ciò mette in evidenza il carattere asintotico 
delle serie che stiamo considerando. Qui l'ordine r di normalizzazione è 
arbitrario, ma lo si può determinare come funzione r = r opt (p) minimizzando 
la funzione (r\) n+1 p r . Si ha così r ~ (1/ ' p) l ^ n+l \ ed utilizzando la formula 
di Stirling si valuta 

|/| ~exp(-(l/p) 1/(n+1) ) . 

Assumendo che il dato iniziale sia contenuto in un intorno A p / 2 dell'origine 
si può allora garantire che l'orbita resterà confinata nel polidisco A p per 
un tempo T ~ exp((l/p) 1 ^ n+1 ^ ) ) che cresce più rapidamente di qualunque 
potenza al decrescere di p. È questo, in forma sintetica, l'argomento che 
conduce alla stima esponenziale del tempo di stabilità tipica della teoria alla 
Nekhoroshev [T6].|17j. 



4 Calcolo effettivo del tempo di stabilità 

Ci poniamo ora l'obiettivo di tradurre l'argomento che abbiamo appena espo- 
sto in un algoritmo di calcolo che ci consenta di dare una valutazione esplicita 
del tempo di stabilità per un sistema reale. Consideriamo un sistema ad n 
gradi di libertà descritto da un'Hamiltoniana della forma ((?]) troncata ad 

2 Uno studio numerico che illustra i meccanismi che conducono alla divergenza si trova 
in [5,0. 
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un ordine r max che può scegliersi, ad esempio, compatibile con la potenza di 
calcolo disponibile. 

Grazie alla manipolazione algebrica costruiamo esplicitamente la forma 
normale di Birkhoff per il nostro sistema fino all'ordine r max . In questa fase 
del calcolo potrebbe presentarsi il problema dei piccoli divisori, ma possiamo 
osservare che grazie al troncamento la condizione di non risonanza deve essere 
verificata solo per \k\ < r max ; questa condizione è facile da verificare dal 
momento che si deve considerare un numero finito di casi. Nel calcolare la 
forma normale avremo cura anche di tener memoria del primo termine del 
resto, ad ogni ordine. In altre parole, per r = 1, . . . ,r max pari costruiamo 
esplicitamente un'Hamiltoniana 

= H (I) +Z 1 (I) + ... + Z<r\l) + F^ r+1 \x,y) + ... + F^\x,y) , 

dove le funzioni F denotano la parte non ancora normalizzata. Di fatto 
la condizione di non risonanza implica che le funzioni Zj si annullino per j 
dispari, ma ciò non ha grande rilevanza per la discussione di questo paragrafo. 

Consideriamo poi un intorno dell'origine a forma di polidisco con raggi 
Ri, ... , R n , ossia 

A pR = {(x,y)eR 2n : x 2 + y 2 <p 2 R 2 , 1 < j < n) . 

Scriviamo un polinomio generico di grado s come 

f(x, y) = J2 fi,kX j y k , \j\ + \k\ = s , 

3,k 

dove abbiamo usato la notazione multiindice j = (ji, . . . ,j n ), k = (ki, . . . , k n ) 
e xiy k = xf ■ . . . ■ y^ k . Scegliamo n parametri positivi R = (Ri, . . . , R n ) e 
calcoliamo la quantità 



In tal modo per p > assegnato abbiamo 

sup \f(x,y)\ < p s \f\ R . 

Questa stima richiede qualche giustificazione. Se consideriamo un disco di 
raggio Ri nel piano X{, yi verifichiamo subito che si ha \xjyf | < R{ +k @j^- 
Basta infatti scrivere Xi = Ri cos 9 , y^ = Ri sin 9 e verificare che sull'intervallo 
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< 9 < 27T si ha | cos- 7 6* sin fe 9\ < Qj,k- La quantità \f\n definita dalla (fTUj) è 
la somma di tutti questi contributi. 

Dobbiamo ora valutare sup^ x ^ eApR \l(x,y)\. Ricordando che la deriva- 
ta temporale di una funzione è la parentesi di Poisson con l'Hamiltoniana 
facciamo uso della diseguaglianza 

\Ì ó {x,y)\<Cp^\{I h F^}\ R , (11) 

avendo scelto una costante C > 1 opportuna. Qui è necessaria qualche 
precisazione perché in linea di principio dovremmo tener conto di una serie 
infinita, il che è chiaramente impossibile in pratica. L'argomento è il seguen- 
te. Dalle stime teoriche sappiamo che la serie dei resti è stimata da una serie 
geometrica. Se p è inferiore al raggio di convergenza della forma normale 
all'ordine r allora esiste una costante C per cui vale la stima riportata sopra. 
Nel calcolo pratico, dal momento che consideriamo p abbastanza piccolo, sce- 
glieremo C = 2. Possiamo però osservare che in pratica la dipendenza del 
risultato finale dalla scelta di C risulta essere molto debole. 

Veniamo dunque al tempo di stabilità. Osservando che Ij < p 2 R 2 /2 
abbiamo anche ìj = R 2 pp, e possiamo riscrivere la diseguaglianza (|TTjl come 

P<^P r+ \ B r , = C\{I v F^}\ R . 

Da qui possiamo ricavare una maggiorazione per la funzione p(t) risolvendo 
l'equazione differenziale p = B r jp r+2 / R 2 . Separando le variabili otteniamo 
che il tempo necessario per passare dal valore iniziale po ad un p > arbitrario 
soddisfa la diseguaglianza \t\ > r(p ,p,r), dove 

. R) ["da . R) ( 1 1 \ 
t Po, Pi t) — min — — / — = min — n- 

Qui introduciamo una scelta per p ponendo p = 2p , per cui la formula 
precedente diventa 

/ 1 \ R 2 
t(pq, 2p , r) = min 1 - — T • 

Poiché questa stima è valida per qualunque scelta di r possiamo ottimizzar- 
la scegliendo il valore r opt (p ) che massimizza r(p ,2p ,r) ed in tal modo 
determiniamo un tempo 

T(po) = max r(p , 2p , r) . (12) 

r 

Questa è la miglior indicazione fornita dal nostro algoritmo, e chiameremo 
T(po) il tempo di stabilità. Osserviamo che tutte le quantità scritte sono 
calcolabili esplicitamente. 
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Tabella 1: Masse ed elementi orbitali eliocentrici di Giove e Saturno calcolati 
dal JPL per il Giorno Giuliano (JD) 2451220.5 . Le lunghezze sono in Unità 
Astronomiche (UÀ); i tempi in anni; la costante gravitazionale è G = 1 . In 
queste unità la massa del Sole è 4n 2 . 





Giove (j = 1) 


Saturno (j = 2) 


massa rrij 
semiasse maggiore a,j 
anomalia media £j 
eccentricità ej 
argomento del perielio Uj 
inclinazione ij 
longitudine del nodo Qj 


(4tt 2 )/1047.355 

5.20092253448245 

6.14053316064644 

0.04814707261917873 

1.18977636117073 

0.006301433258242599 

3.51164756250381 


(4tt 2 )/3498.5 

9.55716977296997 

5.37386251998842 

0.05381979488308911 

5.65165124779163 

0.01552738031933247 

0.370054908914043 



5 Applicazione al sistema Sole— Giove— Saturno 

Veniamo infine all'applicazione al sistema Sole-Giove-Saturno. La prima 
parte del calcolo, del tutto classica, consiste nel calcolare gli sviluppi in 
serie di potenze e trigonometriche necessari per dare una forma esplicita 
all'Hamiltoniana del problema dei tre corpi in variabili di Poincaré, nella 
forma (|2J) e (J3j). Qui scegliamo i valori delle masse e dei parametri orbitali 
di Giove e Saturno, riportati per completezza nella tabella [TJ Nello sviluppo 
teniamo conto dei contributi fino all'ordine 2 nelle masse e fino al grado 6 
nelle variabili lente £, rj. 

Tutto il calcolo è stato svolto grazie ad un pacchetto di manipolazione 
algebrica realizzato ad hoc dagli autori, ed in grado di eseguire tutte le ope- 
razioni algebriche necessarie per l'applicazione degli algoritmi perturbativi. 

Procediamo poi al calcolo della parte secolare del sistema, con un'appros- 
simazione valida fino all'ordine 2 nelle masse, seguendo lo schema illustrato 
nel paragrafo [2J Otteniamo così lo sviluppo dell' Hamiltoniana secolare nella 
forma (j^J) e procediamo al calcolo delle frequenze ed alla diagonalizzazione 
della parte quadratica. Infine calcoliamo la forma normale di Birkhoff fino 
all'ordine 18, che si è rivelato sufficiente per i nostri scopi. 

Avendo calcolato la forma normale di Birkhoff applichiamo l'algoritmo 
di stima del tempo di stabilità descritto nel paragrafo SJ I risultati sono 



illustrati nelle figure l(a) e l(b) Nel primo grafico è riportato il tempo di 
stabilità T(po) in funzione del raggio po nel quale sono contenuti i dati iniziali, 
in scala semilogaritmica. Si nota subito la crescita molto rapida del tempo 
stimato quando p decresce. Che tale crescita sia più rapida di qualunque 
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(b) Ordine della normalizzazione ottimale 



Figura 1: (a) Stima del tempo di stabilità per il sistema Sole-Giove-Saturno 
al variare del raggio po del dominio contenente i dati iniziali. I raggi R per 
il calcolo della norma delle funzioni sono scelti in modo che i dati reali per i 
due pianeti corrispondano a po = 1, e l'unità di tempo è l'anno terrestre. La 
scala verticale riporta il logaritmo decimale del tempo. La linea tratteggiata 
orizzontale indica l'età stimata del sistema solare, (b) L'ordine ottimale di 
normalizzazione r max in funzione di po- 
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potenza lo si arguisce osservando il secondo grafico, in figura 1 (b) in cui 
si riporta l'ordine ottimale r opt in funzione di Pq. I vari tratti orizzontali 
corrispondono ad intervalli in cui r opt (po) resta costante, ed in ciascuno di 
questi tratti il tempo cresce come una potenza p r ° pt - Ma al decrescere di po 
l'ordine ottimale r opt cresce, ed anche molto rapidamente: nel nostro grafico 
il limite r max = 18 viene raggiunto per p — 0.7, e poi resta costante solo 
perché non abbiamo spinto il calcolo ad ordini più elevati. 

Nel grafico l'unità di tempo è l'anno terrestre, ed i raggi R sono scelti in 
modo che i dati correnti per i due pianeti si trovino nel polidisco di raggio 
Po = 1- La retta tratteggiata orizzontale indica l'età stimata del sistema 
solare, corrispondente a circa 5 x IO 9 anni. Dalla figura si vede che il tempo 
di stabilità stimato col nostro metodo risulta essere di circa 1.5 x IO 7 anni. 



6 Commenti e possibili sviluppi 

La stima che abbiamo ottenuto si rivela ancora pessimistica, in contrasto 
ad esempio con le simulazioni numeriche che danno tempi molto più lunghi 
anche per il sistema dei quattro pianeti maggiori. Possiamo però osservare 
che non siamo terribilmente lontani dall'obiettivo, che sarebbe ragionevole, 
di raggiungere almeno l'età del sistema solare. Dal grafico, ad esempio, si 
vede che basterebbe che l'eccentricità fosse pari a poco più di 0.7 volte quella 
reale. La domanda che si pone spontaneamente è se si possano migliorare i 
nostri risultati. 

Una prima osservazione è che lo schema di calcolo che abbiamo seguito 
non è esente da approssimazioni che possono avere un peso rilevante. In effetti 
l'algoritmo di stima dei tempi di stabilità suppone, in buona sostanza, che 
la perturbazione agisca sempre in modo da incrementare l'eccentricità. Ciò 
è certamente falso, ma è difficile tenerne conto nelle stime semianalitiche, 
mentre le simulazioni numeriche svolte mediante integrazione diretta delle 
equazioni del moto ne tengono conto, di fatto. In questa luce, il risultato 
da noi ottenuto può già considerarsi apprezzabile, ed in effetti si colloca tra 
i migliori che vengono tipicamente ottenuti quando si fa ricorso a metodi 
perturbativi. 

Si pone però un problema più profondo, che chiama in causa in modo 
diretto il metodo di Lagrange. Come abbiamo avuto modo di osservare, il 
calcolo delle frequenze secolari è stato svolto da Lagrange facendo riferimen- 
to alle orbite circolari. Se però poniamo il problema della stabilità su tempi 
molto lunghi l'approssimazione dell'orbita circolare può rivelarsi troppo roz- 
za. In effetti sappiamo bene che le orbite circolari non sono soluzioni delle 
equazioni del problema dei tre corpi, e proprio per questo abbiamo cerca- 
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to un'approssimazione migliore che tenesse conto anche delle perturbazioni 
fino al secondo ordine nelle masse. Il nostro calcolo mostra che ciò non è 
sufficiente: le eccentricità del sistema reale sono ancora troppo alte. 

Ci si chiede allora se si possa migliorare il risultato prendendo come rife- 
rimento delle orbite che abbiano già eccentricità lontana dallo zero. Un tal 
procedimento in linea di principio è possibile. In effetti abbiamo già mostrato 
in una memoria precedente [6] che nelle vicinanze dei dati iniziali di Giove e 
Saturno esistono soluzioni quasiper io diche del tipo descritto dal teorema di 
Kolmogorov [8]. Tali soluzioni possono ben esistere anche se si includono nel 
modello i quattro pianeti maggiori, ma dobbiamo ricordare che per sistemi 
a più di 2 gradi di libertà l'esistenza di soluzioni quasiper io diche non è suf- 
ficiente a garantire la stabilità. Siamo quindi portati ad indagare l'esistenza 
di un intorno dell'orbita quasi periodica che sia stabile per tempi molto più 
lunghi di quelli che abbiamo stimato in questa nota. In effetti tale calcolo 
è stato svolto nel lavoro [7], ma si scontra con la difficoltà pratica di dover 
calcolare un numero troppo elevato di termini e non è raggiungibile con la 
potenza dei calcolatori attualmente a nostra disposizione. 

Non resta quindi che cercare approssimazioni migliori dell'orbita circolare, 
ma non eccessivamente impegnative dal punto di vista del calcolo. Questo è 
un problema aperto al quale stiamo dedicando i nostri studi. 
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